Mon. Not. R. Astron. Soc. OQO.rilfIol(2009) Printed 4 January 2010 (MN KT^ style file v2.2) 



Tight constraints on the existence of additional planets around 
HD 189733 ^ f 



O" M. Hrudkova^ ^'^t 1. Skillen^, C. R. Benn^, N. P. Gibson^ ^ D. Pollacco^ D. Nesvorny^ 
O : T. Augusteijn^, S. M. Tulloch^ and Y. C. Joshi^'^ 

, ^Astronomical Institute of the Charles University, V Holesovickdch 2, CZ-I80 00 Praha 8, Czech Republic 
• ^Isaac Newton Group of Telescopes, Apartado de Correos 321, E-387 00 Santa Cruz de la Palma, Canary Islands, Spain 
S ' ^Thiiringer Landessternwarte Tautenburg, Sternwarte 5, D-07778 Tautenburg, Germany 
^ ' * Queen's University Belfast, University Road, BT7 INN, Belfast, UK 
^School of Physics, University of Exeter, EX4 4QL, Exeter, UK 

^Department of Space Studies, Southwest Research Institute, 1050 Walnut Street, Suite 400, Boulder, CO 80302 
^Nordic Optical Telescope, Apartado de Correos 474, E-387 00 Santa Cruz de la Palma, Canary Islands, Spain 
' ' ^Aryabhatta Research Institute of Observational Sciences, Manora Peak, Nainital-263 129, India 

Oh 
6 



4 January 2010 



o 
p 

o 
o 



ABSTRACT 

We report a transit timing study of the transiting exoplanetary system HD 189733. In total 
we observed ten transits in 2006 and 2008 with the 2.6-m Nordic Optical Telescope, and two 
transits in 2007 with the 4.2-m William Herschel Telescope. We used Markov-Chain Monte 
Carlo simulations to derive the system parameters and their uncertainties, and our results 
are in a good agreement with previously published values. We performed two independent 
analyses of transit timing residuals to place upper mass limits on putative perturbing planets. 
The results show no evidence for the presence of planets down to 1 Earth mass near the 1 :2 
and 2:1 resonance orbits, and planets down to 2.2 Earth masses near the 3:5 and 5:3 resonance 
orbits with HD 189733b. These are the strongest hmits to date on the presence of other planets 
in this system. 
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1 INTRODUCTION 

Ground-based radial velocity and photometric transit surveys have 
proved to be the most successful methods for discovering exoplan- 
ets over the past decade, yielding more than 400 extrasolar planets 
discovered to datJB- Most of the exoplanets detected are of Jupiter 
mass, but Earth-mass planets remain to be found. An additional 
planet in a transiting system will perturb the motion of the tran- 
siting planet, and the interval between the mid-eclipses will not be 
constant. Deviations from the predicted mid-transit times can there- 
fore reveal the presence of other bodies in the system, or place lim- 
its on their existence. Short-term variations can uncover the exis- 



* Based on observations made with the Nordic Optical Telescope, operated 
on the island of La Palma jointly by Denmark, Finland, Iceland, Norway, 
and Sweden, in the Spanish Observatorio del Roque de los Muchachos of 
the Instituto de Astroffsica de Canarias. 

f Based on observations made with the William Herschel Telescope oper- 
ated on the island of La Palma by the Isaac Newton Group in the Spanish 
Observatorio del Roque de los Muchachos of the Instituto de Astroffsica de 
Canarias. 

X E-mail: marie@tls-tautenburg.de 

^ The Extrasolar Planets Encyclopedia: |http://exoplanet.eu| 



tence of other planets ('Agol et al."2005VHolman & M urravll2005b . 
moons (Sartoretti & Schneider 1999; Kipping 2009) and also Tro- 
jans Jpord & Holmanj | 2007h. wherea s long-term variations result 
from orbital decay ( Rasio et ai]|l996h and from orbital precession 



induce d by another planet, stel l ar oblateness and genera l relativistic 
effects ('Miralda -Escudil2002l : lHevl & Gladmanll2007t) . Discovery 
of additional bodies can constrain theories of planetary system for- 
mation and evolution. In this paper we describe a transit timing 
study of the transiting exoplanet system HD 189733. 

The HD 189733 transiting system is one of the best studied 
systems from the ground. HD 189733 is a bright star with mag- 
nitude V=l .61 which is orbit ed by a tr ansiting Jupiter-mass planet 
in a period of ~ 2.22 days teouchv et al. 20 05i), and which als o 
has a distant mid-M dwarf binary companion dBakos et al.ll2006al) . 
In 2006 HD 189733 was observed with the M05r (Microvariabil- 
ity and Oscillations of STars) satellite and these data were used 
to search f or the existence of other bodies in the system. First, 
ICroU et al. l (2007) searched for transits from exoplanets other than 
the known hot Jupiter, with the result that any additional close-in 
exoplanets on orbital planes near that of HD 189733b with sizes 
ranging from about 1.7 - 3.5 R®, where is the Earth radius, are 
ruled out. Second, an analysis of transit timing variations (TTVs) 
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Table 1. Observations of the HD 189733 system. The UT date is the date of the beginning of each night. The cycle number is in periods from the ephemeris 
given by^gol et al. 1 2009). For some nights the exposure time was changed during the observations; this is indicated by the second value in parentheses. The 
data rms is per exposure for the ratio of intensities of the target and the comparison star. The barycentric mid-transit times of the HD 1 89733 system are given 
with uncertainties defined as 68 per cent confidence limits. 



Telescope 


UT date 


Cycle no. 


CCD window size 


Exposure 


Data rms 


Mid-transit time 


O-C 


Comment 








(pixels) 


(s) 


(mmag) 


(BJD - 2450000) 


(s) 




NOT 


2006 July 18 


-155 


[1040:200] 


2.5 


2.9 


3935.55805 ± 0.00028 


38 ±25 




NOT 


2006 August 07 


-146 


[1040:200] 


2.5 


2.6 


3955.52509 ±0.00014 


26± 12 




NOT 


2006 August 27 


-137 


[1040:200] 


2.5 (3.0) 


2.7 


3975.49194 ±0.00021 


-2 ± 18 




WHT 


2007 August 17 


-1-23 


[1071,546] 


10.0 


4.6 


4330.46305 ± 0.00042 


-79 ± 36 




WHT 


2007 September 17 


-1-37 


[1071,546] 


3.0(3.5) 


4.4 


4361.52352 ±0.00044 


-43 ± 38 




NOT 


2008 June 07 


-1-156 


[1040:200] 


3.5 (3.0) 


2.6 


4625.53404 ± 0.00038 


-35 ± 33 


partial 


NOT 


2008 June 18 


-1-161 


[1040:200] 


3.5 


2.3 


4636.62768 ± 0.00018 


31 ± 15 




NOT 


2008 July 08 


-1-170 


[1040:200] 


3.5 (4.0) 


2.4 


4656.59451 ±0.00012 


1 ± 11 




NOT 


2008 July 17 


-1-174 


[1040:200] 


3.5 (3.0) 


3.4 


4665.46951 ± 0.00029 


62 ± 25 




NOT 


2008 July 28 


-1-179 


[1040:200] 


3.0 


2.3 


4676.56188 ± 0.00019 


18± 16 




NOT 


2008 August 26 


-1-192 


[1655:200] 


3.5 


2.7 


4705.40332 ± 0.00019 


15 ± 16 




NOT 


2008 September 15 


-1-201 


[1655:200] 


2.5 


2.9 


4725.37064 ± 0.00053 


28 ±46 


partial 



in these data has been carried out bv lJVIiller-Ricci et al] JlOOsI) who 
found that there are no TTVs greater than ±45 s, which rules out 
planets of masses larger than 1 and 4 M^, where is the Earth 
mass, in the 2:3 and 1 :2 inner resonances, respectively, and planets 
greater than 20 in the outer 2: 1 resonance of the known planet 
and greater than 8 IvI^ in the 3:2 resonance. 

Analyses of transit times similar to lMiller-Ricci et al] J2008l) 
have been carried ou t for other transiting planetary systems. 
ISteffen & Agol ( l2005h found no evidence for a second planet in 
the TrES-1 system, excluding planets down to Earth mass near the 
low- order, mean-motion r esonances of the transiting planet. Simi- 
larlv. lGibson et al" Il l2009afbl) found no evidence for additional plan- 
ets down to sub-Earth masses in the interior and exterior 2: 1 reso- 
nances of the TrES-3 and HAT-P-3 systems. 

To measure times of mid-transits with sufficient accuracy to 
detect terrestrial mass planets requires high quality photometry, 
free o f systematic effe cts. HD 189733 is known to have surface 
spots: [Pont et alj feoOTi) observed two spot events in HST (Hubble 
Space Telescope) data when the flux during the transit changed by 
1 and 0.4 mmag. The presence of surface spots on HD 189733 
complicates any transit timing analysis (IVIiller-Ricci et al. 2008). 
The light curve can be distorted if the planet transits in front of a 
spot or due to intrinsic variability of the star. The system parame- 
ters and the mid-eclipse times derived can then be affected by an 
inappropriate fitting model. 

Instrumental effects during transit ingress or egress can also 
influence the accuracy and determination of transit times. For ex- 
ample, if the transit light curve is not properly normalised so that 
all data points in egress have a flux level that is slightly too high, 
the transit time will be determined too early. Correct normalisation 
is especially problematic for partial transit light curves. Both in- 
strumental effects and stellar variability can cause that a light curve 
is improperly normalised. 

It is also important to have a light curve that is well-sampled 
during both ingress and egress, because the transit timing informa- 
tion is contained in these parts. When using large telescopes for 
such a bright star, only short exposure times are needed to get suf- 
ficient signal-to-noise and to avoid saturation, and so the cadence 
of observation is higher. For a given data accuracy, higher cadence 
leads to more accurately determined transit times. 

In section fj2]we describe our observations, and in section ij3] 



we present our data reduction. In section ij4]we explain the tech- 
niques used to estimate uncertainties in our data and to measure the 
system parameters. Finally, in section ij5]we describe the 3-body 
simulations used to place limits on the existence of other bodies in 
the system, and we conclude and discuss our results in section ^ 



2 OBSERVATIONS 

We observed eight full and two partial transits of HD 189733 with 
the 2.6-m Nordic Optical Telescope (NOT), La Palma, Spain, using 
ALFOSC (the Andalucia Faint Object Spectrograph and Camera), 
and two full transits using the AG2 camera on the 4.2-m William 
Herschel Telescope (WHT) of the Isaac Newton Group (ING), La 
Palma, Spain (Table[Tll. 

ALFOSC has a 2048 x 2048 back-illuminated CCD with scale 
0.19 arcsec/pixel and field of view (FOV) 6.5 x 6.5 arcmin^. To 
reduce the readout time of each exposure and the duty cycle of 
observation we windowed the CCD with the window sizes sum- 
marized in Table [T| We used a Stromgren y filter to minimize ef- 
fects of colour-dependent atmospheric extinction on the differen- 
tial photometry and the effect of limb darkening on the transit light 
curves. We defocused the telescope typically to 3.4 arcsec, spread- 
ing the light inside full width at half maximum (FWHM) of the 
Point Spread Function (PSF) over ^ 250 pixels, in order to mini- 
mize the impact of pixel-to-pixel sensitivity variations, and to pre- 
vent saturation. Exposure times were chosen to keep counts below 
50,000 per pixel to avoid saturation of features such as hot spots 
and speckles in the defocused stellar images, and to ensure data 
linearity. The typical exposure time for the NOT data was 3 s (Ta- 
ble[B. 

AG2 is a frame-transfer CCD mounted at the WHT's folded 
Cassegrain focus, based on an ING-designed autoguider head. The 
FOV is 3.3 X 3.3 arcmin^ and the scale is 0.4 arcsec/pixel. We 
used a Kitt Peak R filter and defocused the telescope to 10 and 
12 arcsec for the two nights, spreading the FWHlVI-light over ^ 
490 and 700 pixels, respectively. The corresponding exposure times 
were 10 and 3 s. 

The mid-time of each exposure was converted to the Barycen- 
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trie Julian Date (BJD) using the program BARCOF0. We use BJD 
throughout this paper, because for this system the Heliocentric Ju- 
lian Date would accumulate an error of up to 4 seconds. 



3 DATA REDUCTION 

Bias subtraction, flat-field correction and aperture photometry was 
performed using standard IRAllfl procedures. 

To ensure a signal-to-noise in excess of 1,000 in our 
Stromgren y-filter flat fields for the NOT data we generated a mas- 
ter flat field for each night using individually weighted normalised 
flat fields from the entire observing season combined with weights 
W — 1 — D/ S, where D is the time interval between each night 
and date of observation, and S is the season length. Applying flat- 
field corrections has only a minor effect on the resulting NOT pho- 
tometry, because of the heavily defocused PSF. 

For the WHT data we determined master flat fields with a 
signal-to-noise greater than 1,000 for both nights. However, we 
identified a position-angle dependent scattered light component in 
the flat fields, which introduced systematic noise in our WHT pho- 
tometry. Therefore we did not apply flat-field corrections. 

We used the star 2MASS 20003818-1-2242065 as our com- 
parison star for the WHT data. In our NOT data there are 
two available comparison stars, 2MASS 20003818-1-2242065 and 
2MASS 20003286-1-2241118. We found the ratio of their mea- 
sured intensities varies by a few mmag with time, as the tele- 
scope tracks across the meridian. This variation correlates with 
small drifts in the positions of the stars on the CCD, suggest- 
ing that some light is being lost from the aperture around one of 
the stars due to the wings of the PSF drifting out of that aper- 
ture. A similar variation is seen for the ratio of the intensities of 
2MASS 20003286-1-2241 1 18 and out-of-transit HD 189733, but not 
for 2MASS 20003818-1-2242065 and HD 189733, suggesting that it 
is light from 2MASS 20003286-1-2241 118 which is being lost. This 
star is the farther of the two from HD 189733, and we conclude 
that the variation in measured intensity is due to a combination of 
the small drifts in stellar position on the CCD, and the variation 
of the defocused PSF across the FOV. We therefore used only the 
comparison star which is closer to HD 189733. 

We used circular, equal diameter, photometric apertures for 
both HD 189733 and the comparison star. A range of aperture sizes 
was tried and that producing the minimum noise in the out-of- 
transit data was adopted and fixed during each night. The aperture 
radius for all stars ranged from 18-29 pixels for different NOT 
nights and the typical FWHM was around 18 pixels (3.4 arcsec). 
For the two WHT nights the aperture radius was 28 and 30 pixels, 
respectively, and the corresponding FWHM was 25 and 30 pixels 
(10 and 12 arcsec). 

We ensured the apertures tracked small drifts in the stellar 
positions on each image by using a large centroiding box of size 

4 X FWHM. During each night drifts in the stellar positions on 
the CCD were less than 7 (NOT) and 4 pixels (WHT). The sky 
background was subtracted using an estimate of its brightness de- 
termined within an annulus centred on each star with a width of 
10 pixels. For each night, differential photometry was computed by 



^ http://sin'ah.troja.mff.cuni.czrmary 

The Image Reduction and Analysis Facility (IRAF) is distributed by the 
National Optical Astronomy Observatories, which are operated by the As- 
sociation of Universities for Research in Astronomy, Inc., under cooperative 
agreement with the National Science Foundation. 



taking the ratio of counts from HD 189733 to the counts from the 
comparison star. We normalised our data using linear fits that were 
computed together with other system parameters as described in ij4] 
The normalised unbinned NOT light curves and binned WHT 
light curves, averaged into 10-second bins to have the similar ca- 
dence as the NOT data, are shown in Fig. [T] along with their best- 
fitting models, residuals and data error bars, as derived in 21 



4 LIGHT-CURVE MODELLING 

To estimate the system parameters we used a parametrized model 
where we assumed a circular orbit around the centre of mass to cal- 
culate the normalised separation, z, of the p lanet and star centres as 
a function of time. The analytic formulas o f lMandel& Agol l l2002h 
were used to calculate the fraction of the stellar flux occulted by the 
planet using z and the planet-to-star radius ratio, p. We assumed a 
quadratic limb darkening law: 



1 - 1tl(l - ^t) — U2{1 - /i) 



(1) 



where / is the intensity, fj, is the cosine of the angle between the 
line of sight and the normal to the stellar surface, and ui and U2 are 
the linear and quadratic limb darkening coefficients. For the NOT 
data we allowed the limb darkening coefficients to be free parame- 
ters, in order to include possible errors in the limb darkening coef- 
ficients into our final system parameters and mid- transit times. For 
the WHT data we adopted valu es ui = 0.4970 and U2 = 0.2195 
from the tables o f IClaretl l lioool) and fixed them in the subsequent 
analysis. These correspond to the Johnson R filter which has similar 
characteristics as the Kitt Peak R filter used. 

To compute our model we folded all the NOT light curves of 
full transit except for the night 2008 July 17 which displays obvi- 
ous systematic changes during transit. In our photometry we cannot 
easily distinguish spot effects from systematic instrumental errors; 
to do so would require the instrumental systematic noise to be much 
less than the predicted spot signatures. We fitted simultaneously 
planetary and stellar radius, Rp and _R,t, respectively, the orbital 
inclination, i, two limb darkening coefficients, ui and U2, transit 
time, 7o,n, and additional two parameters for each night n - the 
out-of-transit flux, foot.n, and a time gradient, torad.n- These two 
parameters were allowed to be free to account for any normalisation 
errors in the data. For each change of i?*, the stellar mass, was 

1/3 

recomputed using the scaling relation i?* cx M^, . W e fixed the 
planetary mass value Mp = 1.15 ± 0.04 Mj dBouchv e t al. 20oJL 
adopted a period P = 2.21857503 ± 0.00000037 d dAgol et 
I2009i) . and using Kepler's third law we updated the orbital semi- 
major axis for each choice of M^,. 

W e ran Markov-Chain Monte Carlo (MCMC) simulations 
jTegm ark et al. 2004; Ford 2006; Holman et al. 2006) with the 
Metropolis-Hastings algorithm ( IFordll2005l) to estimate the best- 
fitting parameters and their uncertainties. From an initial point, a 
chain is generated by iterating a jump function, which adds a ran- 
dom value selected from a Gaussian distribution with mean and a 
standard deviation 1, scaled by a factor specific for each param- 
eter so that 44 per cent of e ach parameter sets are accepted 
jGelman et all I2OO3I; iFord l2006h . In each step of the generated 
chain the fitting statistic for old and new parameter values is 
computed: 



NoOF 

E 



fijobs) - fi{theor) 



+ 



(2) 



'Mo 
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Figure 1. Differential piiotometry of tlie HD 189733 system over-plotted with the best-fitting model (solid line) from the MCMC fit. The residuals and Icr 
error bars are also plotted, offset by a constant flux for clarity. The phase was computed using best-fitting transit times presented in Table [T] The photometry 
for NOT data is unbinned and for WHT is binned in time with 10 second bins to give the similar cadence as for NOT data for clarity. 



Limits on additional planets in HD 189733 5 



Here fiipbs) is the flux observed at time i, at is the correspond- 
ing uncertainty, fi(/:/ie or) is the flux calculated using formulas of 
iMandel&Agol (12002') and Ndo f is the number of measurements 
in each light curve. The new parameter is then accepted if its is 
lower than that for the previous parameter, or accepted with a prob- 
ability p = exp ( — if its is higher. The second term in 
Eq. l|2j is a Gaussian prior placed on Ah, where Afo — 0.82 M0 
and = 0.03 Mq is the stellar mas s and its u ncertainty, es- 
timated from stellar spectra by Bouc hv et al.l l l2005h . This ensures 
that errors in the stellar mass, which are the greatest source of un- 
certainty when deriving the system parameters and transit times, 
are taken into account. 

The scale factors were chosen so that ^ 44 per cent o f pa- 
rameter sets were accepted dGelman et alj|2003l : lFordll2006h . For 
each simulation we created 10 independent chains, with length at 
least 100,000 points per chain to ensure convergence. Each chain 
was initiated by a parameter that was within ±5o- of a previously 
known best-fitting parameter value using the estimated uncertainty 
a. The first 20 per cent of each chain was discarded to minimize 
the effect of the initial conditions. We checked c onvergence of gen- 
erated chains using the lOehnan & RubinI l [T992h R statistic and cre- 
ated chains until R < 1.03, a good sign of convergence. 

To estimate appropriate error bars in our data accounting 
for any correlated nois e, we used a proce dure similar to that of 
lOillon et all j2006h and lNarita et alj ( l2007h . We assigned the same 
error bars to all the data points including only Poisson noise. An 
initial MCMC analysis of the folded NOT light curves was used to 
estimate the parameters Rp, i?*, i, ui, U2, To,n, foot,n and tcrad.n- 
The first model light curve was used to find the differences be- 
tween the data and the model for each individual night. Then we 
rescaled the error bars to satisfy the condition /Ndof ~ 1.0, 
where Ndof is the number of measurements in each light curve. 
For the night of 2008 July 17, the nights of the two partial tran- 
sits (2008 June 07 and 2008 September 15) and for the WHT light 
curves (2007 August 17 and 2007 September 17) we adopted our 
first model and ran MCMC analysis to find initial parameters To, 
foot and torad for each night independently. Then we rescaled the 
error bars similarly as before. We assume that our initial model is 
a good description of the light curve. Compared to this model, we 
found that for the NOT data errors are higher by factors of 2.3 - 3.4 
than errors including only Poisson noise, and for the WHT data by 
factors of 4.6 and 4.4 for the two nights, respectively. The data rms 
errors per exposure are presented in Table[T] The predicted rms due 
to photon noise, which is dominated by the fainter comparison star, 
and to atmospheric scintillation, is ~ 2.5 mmag for the NOT data 
and ~ 3 mmag for the WHT data. 

The amplitude of systematic trends in the photometry was es- 
timated from the standard deviation over one residual point, (Ji, 
and from the standard deviation of the average of the residuals over 
A*' successive poin ts, ajy. We solved the following system of two 
equations given bv lGillonetal] ( l2006t) : 

fl^fm+cr?, (3) 

2 

0-iV = ^ + (^r, (4) 

to obtain the amplitude of the white noise, ct„, which is uncorre- 
lated and averages down as and the red noise, (Jr, which 

is correlated and remains constant for specified A^. The error bars 
were then adjusted by multiplying by [1 + 7V(cr, /(7u,)^]^'''^ and 
these rescaled uncertainties were used for the subsequent fitting 
procedure. To account properly for the systematic errors, the result- 
ing multiplying factor was computed as the average of values using 



Table 2. System parameters of HD 189733. The uncertainties are 68 per 
cent confidence limits. 



Parameter 


Symbol 


Value 


Units 


Planet radius 


Rp 


1.142 ±0.014 


Rj 


Star radius 


R* 


0.755 ±0.009 


Rq 


Orbital inclination 


i 


85.70 ±0.03 


deg 


Planet/star radius ratio 


P 


0.1556 ±0.0027 




Total transit duration 




1.807 ±0.023 


h 


Impact parameter 


b 


0.667 ± 0.009 





different in the range 15 - 30 minutes (the typical time-scale of 
ingress and egress). 

To create our final model, we proceeded as before but this time 
including systematic noise in our data and therefore properly esti- 
mating parameter uncertainties. We ran MCMC using the folded 
NOT light curves and fitting the parameters as described earlier. 
We created 10 chains, each with length 2,000,000 points in order 
to achieve convergence. Ultimately, we used our final model to find 
individual mid-eclipse times and two normalisation parameters for 
the night of 2008 July 17, the nights of the two partial transits (2008 
June 07 and 2008 September 15) and for the WHT light curves 
(2007 August 17 and 2007 September 17). 



5 RESULTS 

The final system parameters are presented in Table |2] and are 
consiste nt within ~ 2a err o r bars with the p r eviously published 
values (iBakose tal. 2006b|; iPont et al.l 120071; IWinn et alj |2007| : 
iMiller-Ricci et al. 2008) . The resulting limb darkening coefficients 
for the NOT data were ui = 0.46 ± 0.10 and U2 = 0.35 ± 0.13. 

The final barycentric transit times can be found in Table [T] 
The uncertainties are defined as 68 per cent confidence limits. To 
compute the observed -minus-calculated values (O — C) we used 
the ephemeris given bv lAgol etalj ( l2009l) : 

Tc{E) = HJD (2454279.436741 ± 0.000023) + (5) 
(2'!21857503 ± 0'?00000037) x E. 

The resulting O — C residuals together with all the other previously 
publi s hed values jBakos et a l.ll2006bl :IPont et al.ll2007l: Iwinn et al.l 
I2OO7I : iMiUer-Ricci et alj I2OO8: Knut son et alj l2009h are plotted 
in Fig. [21 Our observations did not bring any refinement of the 
ephemeris and we confirm that presented by A gol et al.„ ^200^. 
For the night 2006 August 07 a transit timing measurem ent of 
HD 189733 was also presented bv lMiller-Ricci et al.l feOOSh from 
the MOST data and it is consistent within 2a error bars with our 
measurement. 

5.1 Transit timing variations analysis 

For all our observations which span more than two years, the mean 
0~C= 5 ± 38 s, where the quoted error is the rms scatter in 
the 0~C values and is slightly larger than the average O — C un- 
certainty ~ 25 s. None of our O — C measurements is a significant 
outlier. The two largest O — C values for the nights 2007 August 17 
and 2008 July 17 coincide with obvious systematic changes during 
the transit (see Fig.[Tll and both have the same or larger uncertainty 
than the average value. Therefore the rms scatter in the O — C val- 
ues of 38 s is a good estimate for placing limits on the presence of 
other planets in the system. 
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Figure 2. Top: NOT and WHT O — C residuals of mid-transit times of tlie HD 189733 system including b oth partial (star symb ol) and full transits (filled circle 
symbol). M iddle: Previously publishe d values p lotted together with NOT and WHT results. Filled squares: Bakos etaljj2006bl) . ground -based; filled tr i angles : 
IPont et alj fcoO?), HST; open squares: lTO"nn et al. ( 2007), ground-based; open circles: Miller-Ricci et al. 1 2008), MOST; open triangles: lKnutson et al.U2009l) , 
Spitzer; filled circles and stars: this work. Bottom: The same as the middle but zoomed for clarity. The cycle number is in periods from the ephemeris given 
by|Agol et al. (2009). A horizontal line is plotted in each panel at O — C = to guide the eye. Our timing measurements are the most accurate from known 
ground-based observations. 
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We used this conclusion to place mass limits on the exis- 
tence of planets on orbits interior and exterior to HD 189733b. 
First, we selected the mass, semimajor axis and eccentricity of 
the putative perturbing planet. The orbital inclination was set so 
that HD 189733b and the perturbing planet have coplanar orbits. 
The two-planet system wa s then numerically integrated using the 
Bulirsch-Stoer integrator jPress etal.lll992l) . We determined all 
mid-transit times of HD 189733b over a time-span of 500 days, 
an interval long enough to cover at least 14 orbits of all perturbing 
planets we can exclude, and used these data to estimate TTVs. The 
mass, initial semimajor axis and initial eccentricity of the perturb- 
ing planet were varied to determine the TTV amplitude for different 
planetary configurations. 

Fig. [3] shows the range of the inner and outer planet's orbits 
that produce TTVs smaller/larger than ±38 s and are thus compati- 
ble/incompatible with our TTV observations of the HD 1 89733 sys- 
tem. The shaded area in Fig.|3]excludes a range of possible eccen- 
tricities and semimajor axes for a putative 1 Earth-mass (top) and 

2 Earth-masses (bottom) inner (left) and outer (right) planet in the 
system. Based on this analysis, our observations of the HD 189733 
system show no evidence for the presence of planets down to 1 
Earth mass in the 2:1, 3:2 and 5:3 exterior resonance orbits, planets 
down to 1 Earth mass in the 1:2, 1:3, 2:3 and 2:5 interior reso- 
nance orbits, and planets down to 2 Earth masses in the 1 :4 interior 
resonance orbit with HD 189733b. However, not all of these res- 
onant orbits a re Hill stable. W e computed Hill stability according 
to Eq. (21) of lGladmanI J 19931) for both inner and outer perturbing 
planet and displayed the result in Fig. |3] using thin solid line. For 
the inner/outer perturbing planet all the orbits on the left/right to the 
thin solid line are Hill-stable, which means that close approaches 
between two planets are forbidden. For the rest of the parameter 
space the Hill stability of the system is unknown; the system still 
may be Hill-stable. 

iNesvorn^l ( l2009 h showed that the TTV signal can be signifi- 
cantly amplified for planetary systems with substantial orbital in- 
clinations of the transiting and perturbing planet and/or in the case 
of transiting planet in an eccentric orbit with an anti-aligned orbit 
of the perturbing planetary companion. Therefore for most orbital 
architectures of exoplanetary systems we determine the perturber's 
upper mass from our TTVs under the assumption of coplanar orbits 
of transiting and perturbing planets. 

However, the above mentioned analysis does not take into ac- 
count time sampling of our measured transit times and their un- 
certainties. It is possible to have a system whose TTV amplitude 
exceeds 38 s but remains consistent with the available transit tim- 
ing data. To assure that the limits on additional planets presented in 
this paper are not overestimated, in addition to our previous analy- 
sis we compared model timing residuals against the transit times to 
place upper mass limi ts for a putative pertu rbing planet. We used 
the same procedure as lGibsonetalJ ( l2009allbl) , where more details 
can be found. To compute model timing residuals we integrated 
the equations of motion for a three body system using a 4th-order 
Runge-Kutta method, with the first two bodies representing the star 
and planet of the HD 189733 system, and the third body represent- 
ing a putative perturbing planet. The transit times were extracted 
when the star and transiting planet were aligned along the direction 
of observation, and the residuals from a linear fit were taken to be 
the model timing residuals. For each model, TTVs were extracted 
for 6 equally spaced directions of observations, and we simulated 

3 years of TTVs to cover the full range of observations. The re- 
sulting TTVs were then compared to transit times presented in the 
middle panel of Fig.|2l i. e. all available transit times. Due to com- 



putational limitations we assumed that the amplitude o f the timing 
residuals scales proportionally to the perturber mass ( lAgol et al.l 
l2005l : lHolman & Murravll2005ll . that a perturber has circular orbit 
and that the orbits of the planets are coplanar. 

We created models with period ratios in the range 0.2 - 5.0, in- 
creasing the sampling around the interior 1 :2 and exterior 2: 1 reso- 
nanc es. The maximum all owed mass for each model was calculated 
as in lGibson et alj j200 9a' b). We scaled the perturb er mass until the 
of the model fit increased by a value Ax^ ~ 9 dSteffen & Agoll 
l2005h from that of a linear ephemeris. Then we minimized the 
along the epoch, and rescaled the perturber mass again until the 
maximum allowed mass was determined. This was repeated for 
each direction of the observation, and the maximum mass found 
was set as our upper mass limit. This process was repeated twice to 
verify our assumption that the timing residuals scale proportionally 
with the mass of the perturbing planet. 

The resulting upper mass limits are plotted as a function of the 
period ratio in Fig.|4] The solid line represents the upper mass limits 
from our three-body simulations, and the horizontal dashed line 
shows an Earth-mass planet. Based on this analysis, the available 
data were sufficiently sensitive to probe for masses as small as 0.2 
and 0.15 M9 near the interior 1 :2 and exterior 2: 1 resonances with 
HD 189733b, respectively. The corresponding upper masses near 
the 3:5 and 5:3 resonances with HD 189733b are 2.2 and 0.54 Mg,. 
In the rest of the space outside the region between the 2:3 and 3:2 
resonances with HD 189733b the upper mass limits are of the order 
of a few tens of Earth masses to a few Jupiter masses. However, 
these upper mass limits result from the assumption of a circular 
orbit of a perturber. Eccentric orbits may lead to smaller TTVs, and 
hence planets larger than our upper mass limits in eccentric orbits 
could exist in these regions. Unfortunately, accounting for eccentric 
orbits is computationally unfeasible using these models due to a 
large parameter space. Thus real upper mass limits of a perturber 
in a low eccentric orbit can be as much as an order of magnitude 
larger ( Gibson et al. 2009b). 

We also consider the possible prese nce of Trojans in the sys- 
tem. According to .Ford & HolmanI l l2007h transit times are the same 
for a system without a Trojan and for a system where the tran- 
siting planet and Trojan have equal eccentricities and the Trojan 
resides exactly at the Lagrange L4/L5 fixed point. TTV analysis 
alone is not suitable for constraining the presence of Trojans in 
transiting systems. However, a comparison of the photometrically 
observed transit time and the transit time calculated from the ra- 
dial velocity data assuming zero Trojan mass can reveal a Trojan 
or place upper limits on its m ass. Such an analysis was done by 
lMadhusu(ttan & Wind ( |2009|) who found an upper limi t of 22 Mm 
for a Trojan in the HD 189733 system. In addition, ICroll et aH 
searched for Trojan transits in MOST photometry, assum- 
ing similar inclinations of the Trojan's and transiting planet's orbits, 
and concluded that Trojans with a radius above 2.7 should have 
been detected with 95 per cent confidence. Using a mean density of 
p ~ 3 000 kg m^^, this corresponds to 11 M®. We used Eq. (1) of 
iFord & Holman (2007) to estimate what Trojan's mass can be ex- 
cluded in the system based on 38 s rms of our TTVs. However, the 
amplitude of the angular displacement of a putative Trojan from the 
Langrange point is not known. If these libration amplitudes are sim- 
ilar as for Trojans orbiting near the Sun-Jupiter Langrange points, 
that is 5 - 30 deg (Murrav & Dermott 2000), our TTVs show no 
evidence for Trojans with masses higher than 5.3 M9. 

For an Earth-mass ex omoon in a circular orbit about 
HD 189733b iKippiiiS ( I2OO9I) predicted TTV amplitude of 1.51 s 
and transit duration variation (TDV) amplitude of 2.94 s. Increas- 
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Figure 3. A numerical survey of the HD 189733 system showing 38-s TTVs caused by an inner (left plots) and outer (right plots) 1 Earth-mass planet 
(m2 = 3 X 10"" M0, top) and 2 Earth-masses planet (m2 = 6 X 10"'^ M©, bottom). The shaded area excludes a range of possible eccentricities and 
semimajor axes for a putative 1 and 2 Earth-masses inner/outer planet in the system based on our observational non-detection of TTVs greater than ±38 s. 
We do not display plots for Jupiter-mass planets as these would easily be detected in radial velocity searches. The thick solid line shows a boundary where a 
collision between the two planets can occur. It is defined so that the apocentre/pericentr e of the inner/out er perturbing planet equals to the semi-major axis of 
the transiting planet. The thin solid line represents a Hill stability computed according to lGladmanI il993h . On the top of the upper panel we indicate the major 
resonances of the putative perturbing planet and HD 189733b. 
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Period Ratio 

Figure 4. Upper mass limits on a putative second planet in the HD 189733 system as a function of period ratio based on the comparison of model timing 
residuals and all available transit times. The solid line represents the upper mass found using three-body simulations. The horizontal dashed line shows an 
Earth-mass planet. The grey area is the region where an Earth-mass perturbing planet is not guaranteed to be Hill stable. 



ing the eccentricity of the moon's orbit decreases TTV amplitude, 
but increases TDV amplitude. However, for the HD 189733 system 
these variations are too small to be detectable in our data. 



6 CONCLUSIONS AND DISCUSSION 

iMiUer-Ricci et alj 12008') found no TTVs greater than ±45 s in 
MOST data, and excluded super-Earths of masses larger than 1 and 
4 in the 2:3 and 1:2 inner resonance, respectively, and planets 
greater than 20 in the outer 2: 1 resonance o f the known planet 
and g reater than 8 in the 3:2 resonance. iMiller-Ricci et al.l 
( 12008 ) assumed that the orbit of the perturbing planet is circular and 
that additional planets i n eccentric orbit s would produce stronger 
perturbations. However. iNesvorti ^l l2009l) showed that an eccentric 
planet can produce stronger or weaker perturbations depending on 
the relative angular position of its orbital pericentre. 

In this paper we used two different methods to determine the 
upper mass limits for a putative perturbing planet in the HD 189733 
system and thus the results of both analyses can be directly com- 
pared. Our first analysis does not take into account time sampling 
of the measured transit times and their uncertainties. On the other 
hand, it was possible to probe for eccentric orbits of a perturb- 
ing planet, whic h is more rigorous than assuming a circular orbit 
In esvomvll2009h . Further analysis was performed to assure that the 
limits on additional planets presented in this paper are not overesti- 
mated. Unfortunately, applying this mothod for eccentric orbits of 
a perturbing planet would increase the number of parameters enor- 
mously, thus we assumed a circular orbit for the perturber. 

Due to the limitations of our TTV analyses, we adopt the 
least constaining limits to conclude what upper masses of a pu- 



tative perturbing planet can be excluded in the HD 189733 system. 
The results show no evidence for the presence of planets down to 
1 Earth mass near the 1:2 and 2:1 resonance orbits, and planets 
down to 2.2 Earth masses near the 3:5 and 5:3 resonance orbits with 
HD 189733b. These are the strongest limits to date on the presence 
of other planets in this system, based on results of two independent 
TTVs analyses. We also discuss the possible presence of Trojans in 
the system, and conclude that the highest limit on a Trojan mass is 
5.3 M0 if its libration amplitude is similar as for Trojans orbiting 
near the Sun- Jupiter Lagrange points. 
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